LARGE SETS OF CONSECUTIVE MAASS FORMS AND 
FLUCTUATIONS IN THE WEYL REMAINDER 
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Abstract. We explore an algorithm which systematically finds all discrete eigenvalues of 
an analytic eigenvalue problem. The algorithm is more simple and elementary as could be 
expected before. It consists of Hejhal's identity, linearisation, and Turing bounds. Using 
the algorithm, we compute more than one hundredsixty thousand consecutive eigenvalues 
of the Laplacian on the modular surface, and investigate the asymptotic and statistic 
properties of the fluctuations in the Weyl remainder. We summarize the findings in two 
conjectures. One is on the maximum size of the Weyl remainder, and the other is on the 
distribution of a suitably scaled version of the Weyl remainder. 



1. Introduction 

The real-analytic eigenvalue problem of the Laplacian on finite, non-compact hyperbolic 
surfaces is an important one. The solutions are automorphic forms and play a central role in 
analytic number theory. They can be used to extend the classical theory of Dirichlet series 
with Euler products [21 J and are closely related to the Millennium Problem of the Riemann 
Hypothesis [29]. Moreover, they serve to find class numbers of real quadratic number fields 
and hence solve questions which already inspired Gauss [5] . Non-holomorphic automorphic 
forms, the so called Maass forms, are intensively studied in spectral theory. The spectral 
decomposition into Maass forms let to the discovery of Selberg's trace formula f30"llll| which 
connects the spectrum to the geometric properties of the underlying space. Automorphic 
forms are the prime example of quantum unique ergodicity [20]. In addition, Maass forms 
span a Hilbert basis in quantum mechanics on hyperbolic surfaces and serve as important 
examples in quantum chaos. For instance, they play a distinguished role in eigenvalue 
statistics [3] . As a complete set of eigenfunctions of the Laplacian on hyperbolic surfaces 
they have also found applications in general relativity and cosmology pp. 

By approximating them, numerics can bring Maass forms closer to us [16]. It is clear 
how to compute Maass forms [13]. However, it was difficult to find consecutive sets of 
solutions. We present an algorithm which allows us to find large consecutive sets of Maass 
forms efficiently. The algorithm is based on three ingredients: Hejhal's identity to compute 
eigenfunctions corresponding to given eigenvalues; linearisation which converts the analytic 
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eigenvalue problem locally to a matrix eigenvalue problem; and explicit Turing bounds 
which serve to check and complete the results. 

As they lie on our route, we apply results of ergodic theory, namely the equidistribution of 
long closed horocycles [21 EH Q21 Ell 033 E2] • This allows us to unreveil an elegant view on 
the derivation of Hejhal's identity in section [4j where we present the computation of Maass 
forms under the preliminary assumption that the discrete eigenvalues would be known. 
The eigenvalue problem is linearised in section [5j thereby almost all of the eigenvalues are 
found. Section [6] presents Turing bounds. These bounds drive a control algorithm which 
systematically checks and completes the list of Maass forms until a large set of consecutive 
Maass forms is found. We have computed more than 160 000 consecutive Maass forms. 
This exceeds any previous numerical solution of any non-integrable system by far. Results 
are listed in section [7J where we investigate the asymptotic and statistic properties of the 
fluctuations in the Weyl remainder. 

2. Fundamental domains, equidistribution of long closed horocycles, and 

congruent points 

In this section, we will consider fundamental domains, present a pullback algorithm, and 
use ergodic properties to show that a discrete subgroup of the isometries can be completely 
specified by a set of congruent pairs of points. 

Let X = T\H be a surface, where T is a cofinite, non-cocompact subgroup of SL(2, K) 
which acts properly discontinuous on the Poincare upper half-plane %. The action is given 
by linear fractional transformations. On X we have the invariant metric ds = \dz\/y. We 
use £(■) to denote length, and d(-, •) to denote hyperbolic distance. 

Because T is non-cocompact, X has at least one cusp. After a suitable coordinate 
transformation we may assume that one of the cusps lies at oo, and that the corresponding 
isotropy subgroup C T is generated by z h- >• z + 1 . 

Take some point pinU, not an elliptic fixed point. Let J- be the corresponding Dirichlet 
fundamental domain: 

T = {z G U | d{p, z) < d{p, 72) V7 G T}. 

We know from Siegel's theorem that J- has a finite number of sides. The sides of J- fall into 
congruent pairs, let us denote them s±, s n+ i; S2, s n+ 2] ■ ■ ■ ; s n , S2 n , and let 51, 52, ■ • • , 9n G T 
be the identification maps, gtSk = Sn+fc- (Recall the convention that an elliptic fixed point 
of order two is considered as a vertex separating two distinct sides of J-.) We know that 
gi,g<z,...,g n generate T (see [TBI PP- 73-74]). 

The following is an algorithm for computing the pullback of any given point in the 
hyperbolic plane into a Dirichlet fundamental domain. 

Algorithm 1 (Pullback algorithm [31] )• (1) Start with z G %. 

(2) Compute the 2n points g\z, g2Z, g^z, . . . , g~ 1 z. Let z' be the one of these 
which has smallest hyperbolic distance to p. 

(3) If d(p, z') < d(p,z), then replace z by z' , and repeat from^ 
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(4) In the other case, d(p,z') > d(p,z), we know that z lies in T, hence z is the 
searched-for pullback. 

Strombergsson proved that his pullback algorithm always finds the pullback within a 
finite number of operations [31j . 

We use z* = x* + iy* to denote the pullback of z = x + iy, and A* = {z* \ z G A} to 
denote the pullback of A, for any A C %. 

For any y > 0, the curve L y = {x + iy \ x G [0, 1]} is a closed horocycle of length t = 1/y 
on X. When t — > oo, this curve is known to become uniformly equidistributed on X with 
respect to the Poincare area d\i = dxdy/y 2 |32| . Let us consider Q equidistant points 
{zj}i<j<Q on the closed horocyle L y . Using AlgorithmJTJ we can compute {z*}\<j<Q, and 
construct the congruent pairs of points {(zj, -zp}i<j<Q- 

Of importance is: 

Lemma 2. Fix Q G N and y > 0, and consider Q equidistant points {zj}\<j<Q on the 
closed horocycle L y . For Q sufficiently large and y sufficiently small, the congruent pairs 
of points {(zj, Zj)}i<j<Q contain all information of the group T. 

For A C Ti, we use dA to denote its boundary, and TA = {^z \ z G A, 7 £ T} to denote 
its orbit. 

For the proof of Lemma [2j we need the following: 

Lemma 3. // the horocycle L* comes sufficiently close to each point in T , then we have: 
For each 1 < k < n there exists a line segment 1^ C L y such that meets TdJ- exactly 
once, l* k does not meet dF in a vertex, and £(h) > £(h H T^T 1 -^ 7 ) > for some 7^. G T. 

Proof. Fix 1 < k < n. Fix some point uj^ G Sk, not a vertex of J 7 . By assumption L* comes 
sufficiently close to any point in T . Thus, for some e > there exists 7^ G T such that 
\y — 9(7 A r 1 w/ c )| < e. If unexpectedly 9(7 s T 1 cjfc) = sup wgSfe ^s(^ 1 oj) should hold, we repeat 
the proof with another ui^ G Sfc . 

We have inf w£Sfc ^{^uj) < Q^ 1 ^^ < y+e and y-e < Qi^Uk) < sup wgSfc 9(7 fc x o;). 
For e sufficiently small, intersects J^s^ at least once (and at most twice) not at a vertex. 

Take an open line segment Ik C L y such that Ik intersects 7^ 1 Sfc, and such that it meets 
TdJ- exactly once, but not in a vertex of 'y^ 1 J-. Since Ik is an open line segment, it has 
positive length on both sides of the intersection with 7^T 1 Sfc. □ 

Proof of Lemma [J| Our task is to extract the side identification maps g\ , . . . , g n from the 
pairs of points {{zj, z*)}i<j< Q . 

The points z\, . . . , zq are on the curve L y . In the limit y — > 0, the closed horocycle L y 
equidistributes on T\H. If y > is finite, the curve L* is no longer dense in J '. Still, for y 
sufficiently small, L* comes sufficiently close to each point in J '. 

Fix some 1 < k < n. Let Ik be the line segment given by Lemma|3j There are two distinct 
elements jk,"fn+k £ T, 7t / ±7n+fc, such that i(l k n 7^ 1 J r ) > and l(l k n > °- 

For Q sufficiently large, at least three succesive points of {^}i<j<Q are on Ik n 7^ 1 -? r and 
another three are on l k Pi J~^ k J-- This is true for any 1 < A; < n. Hence, for each 1 < k < n 
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there are six succesive pairs of points in {(zj, z*)}\<j<Q such that three of them satisfy 
z* = jk z j, while the other three satisfy z* = ^ n+ kZj. 

Now, let only the congruent pairs of points {(zj, ^)}i<j<q be given. We search for all 
sequences of six succesive pairs of points in {(Zj, z*)}i<j<q subject to the condition that 
three of them satisfy z^ = ^hZj, while the other three of them satisfy z*- = "fh'Zj, where 
Ihilh' ^ SL(2, R), 7/j / ±7h'- By construction all 7^ and 7^/ are in T. Moreover, amongst 
all these sequences of six succesive pairs of points, for each 1 < k < n there is at least one 
sequence which is associated to the congruent pair of sides Sk,s n+ k- The corresponding 
side identification map (or its inverse) follows from = Jhjh' ■ '— ' 

3. MAASS FORMS 

Let X = T\H be a finite, non-compact surface with invariant metric ds = \dz\jy. 
According to the metric, the Laplacian reads A = —y 2 (d 2 /dx 2 + d 2 /dy 2 ). 

Definition 4 ([21]). A Maass form is a 

(1) real analytic, f G C°°(H), 

(2) square integrable, f £ L 2 (X), 

(3) automorphic, f{"fz) = f(z) V7 E T, 

(4) eigenfunction of the Laplacian, Af(z) = \f(z). 

If a non-constant Maass form vanishes in all cusps of X, it is called a Maass cusp form. 

According to the Roelcke-Selberg spectral resolution of the Laplacian [301 [2^, its spec- 
trum contains both, a discrete and a continuous part. The discrete part of the spectrum 
is spanned by the constant eigenfunction /o and a countable number of Maass cusp forms 
/i) /2j /3, • • • which we take to be ordered with increasing eigenvalues = Ao < Ai < A2 < 
A3 < The continuous part of the spectrum A > 1/4, is spanned by Eisenstein series. 

In order to keep notation simple, we focus on Maass cusp forms on finite surfaces with 
exactly one cusp. We assume that the cusps lies at 00, and that the corresponding isotropy 
subgroup Too C T is generated by z i->- z + 1. 

In this setting, Maass cusp forms associated to the eigenvalue A = r 2 + 1/4 have the 
Fourier expansion 

/(*)= Yl any l/2 K ir (2K\n\y)e 2 ™ nx , 

nSZ-{0} 

where K stands for the JC-Bessel function which decays exponentially for large arguments 
Eq. (14)]. The expansion coefficients a n grow at most polynomially in n. If we bound 
y from below, y > Y, and allow for a small numerical error of size e > 0, there is an 
M(Y) = const I Y such that 

(1) /(*)= Yj any l/2 K tr (2it\n\y)e 2mnx + [\error\ < e]. 

0¥=H<M(Y) 

The constant in M(Y) = const jY , which depends on e > 0, r, and the group T, can be 
worked out explicitly in each case. Of importance is: 
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Lemma 5. Let the fundamental domain T be bounded from below, Yq := inf zg jr ^s(z) > 0. 
Apart from a numerical error of at most e > 0, a Maass cusp form is completely specified 
by its eigenvalue and a finite number of expansion coefficients, {A, {a n }}o^| n |<M(Y ) - 

Proof. Let {A, {a n }}o^\n\<M(Y ) be given. By automorphy it is enough to know f(z) inside 
the fundamental domain. There, we have y > Yq, and the value of f(z) follows from 
Q. □ 

Remark 6. For any {A, {a>n}}o^\n\<M(Y )> the Fourier expansion ([!]) is real analytic, square 
integrable on X , vanishes in the cusp at oo, and satisfies the eigenvalue equation of the 
Laplacian. However, only for very specific choices of {A, {an}}o^\n\<M(Y )> the correspond- 
ing function is automorphic, and hence a Maass cusp form. 

4. Computing Maass forms 

For the moment, we assume that the eigenvalue would be given. Our task is to find 
the coefficients {a„} such that automorphy holds, f(jz) = f(z) V7 G T. We use the 
congruent pairs of points of Lemma[2]with < y <Yq and Q > 2M(y) to test automorphy, 
f(Zj) = f{zj) for all 1 < j < Q. This results in Q linear equations in 2M unknowns 

{a n }o^\n\<M ■ 

Note, however, that f(zj) = f(z*) is ill-conditioned. If one tries to solve for aM, one 
gets a M y 1 / 2 K ir (2TT\M\y)e 27TlMx ^ = f{z*) - [f{zj)i ess the M _ th summand ]. The absolute value 
of the l.h.s is bounded by e which is smaller than the numerical error of the r.h.s.. 

We use exponentially weighted superpositions to convert the ill-conditioned linear system 
of equations into a well-conditioned linear system of equations |13j . 

h Q -. Q 

Using geometric series on the l.h.s. and Fourier expanding the r.h.s. results in Hejhal's 
identity: 



1 Q 

(2) a m y l l 2 K ir {2^\m\y) + [\error\ < e) = - £ /(^e^"" ' < 



3=1 
1 Q 

-^/(4)e- 2 ™^ = Y, a n V mn + [\error\<e\ V0 ^ |m| < M(y) 

3=1 <tf\n\<M(Y ) 



where V mn = A £? =1 y* 1/2 K ir (2ir\n\y*)e 



2iri(nx*—mxj) 



Algorithm 7 (Phase 2 p2]). //{A, {a n }}o^\n\<M(Y ) are given, then any further coefficients 
follow directly from 

a ™ ~ V i/2 K . \ 2 «\m\y) S V ™ a " V|m|>M(y ) 
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as y — > 0. The numerical error is bounded by 

t i 2e 
\error of a m \ < , , — ; — — — . 

We use C(A) to denote the matrix C mn (A) = ^n- yl/2 ^ 2 ; |m|;;) , ^ \m\ < M(Y ), ^ 
|n| < M(y ), and a to denote the vector a = (a_ M (y ), a_ M (y )+i, • • • , a-i, ai,a 2 , . . . ,a M ( Yo ))- 

Algorithm 8 (Phase 1 [13]). //only A is given, then the coefficients {an}o^|n|<M(y ) f°^ ow 
from solving the linear system of equations 

C(A)a~0, 

The numerical error is bounded by 

error of V C mn a n < . -| V0 ^ \m\ < M(Y ). 

^ \y 1 ' z K ir (2ir\m\y)\ 



0^\n\<M(Y ) 



Remark 9. // we make a good choice for the height y of the horocycle and take some 
Mq > M(Yq) such that \y l / 2 Ki r (2it\m\y)\ 3> 1 for all ^ \m\ < Mq, the linear system of 
equations Q is well-conditioned for 7^ \m\ < Mq. 

5. Finding eigenvalues 

It remains to find the eigenvalues A for which the linear system C{X)a ~ has non-trivial 
solutions a / 0, and for which the non-trivial solutions are independent of the height y of 
the horocycle. 

To this end, we discretise the A-axis into trial values < Ai < A2 < A3 < . . . which lie 
sufficiently close together, and explore the neighbourhood of each trial value for non-trivial 
solutions of C(A + h)a ~ 0. For the latter, we linearise with respect to the eigenvalue, cf. 

h 2 - 

C(A + h) = C(A) + hC'(X) + —R{\, h), 

where R(X, h) is a matrix whose norm is bounded by 

\\R(X,h)\\ < sup \\C"(X + h)\\. 

\t\<W 

Let A be fixed, n be a small complex parameter, and let (h v ,d v )^ be the non-trivial 
solutions of 

2 

(3) C'-\X){C(X) + ^RCX,v))a v = -h v a v . 

For 7] = 0, we know (h, a)a : = Qiq, ao)^ numerically. If C"(A) is regular, there are 2M(Yq) 

matrix eigenvalues of C /_1 (A)C(A), and there are 2M(Yq) generalized eigenvectors. 
For r] 7^ 0, we use perturbation theory to find (h v , a v )^. 
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Conjecture 10. Under suitable conditions, in particular \r]\ sufficiently small, we have 

h ri = h + 0(rj 2 ) and a v = a + 0(ri). 



Apart from numerical confirmations, it would be nice to have a proof of Conjecture 10 
However, it is not even evident that all solutions of ^ have a Taylor expansion in r\ around 
0. 



If Conjecture 10 is true, it implies the following: Using ?j = /i^ in Q, we get C(A + 
r\)a^ = with r\ = h + 0(n 2 ) and oL = a + 0(rj). For \n\ sufficiently small this can be 
solved for n which gives r\ = h + 0(h 2 ). This allows to find numerical approximations, 
A = A + h + 0(h 2 ), for the eigenvalues of C(A)a ~ in some neighbourhood of the trial 
value A. The approximations can be refined iteratively. This establishes: 



Theorem 11. If Conjecture 10 holds true, and if the trial values < Ai < A2 < • ■ • < 

Amax He sufficiently close together, the linearisation in the eigenvalue allows to find all 
eigenvalues ofC(X)a ~ in the interval [0, X max \. 

Remark 12. (1) Not each solution of (C"~ 1 (A)C(A) + /i/)a = approximates a Maass 
form. Clearly, \h\ can be too large such that the linearisation becomes inappropri- 
ate. This is of no concern, since we do not need to consider Maass forms whose 
eigenvalues are further away than the next trial value A. 

(2) It can happen that a formal solution is dependent on y. This has to be checked 
explicitly by re-evaluating C(X)a for different heights y of the horocycle. Only if 
C(X)a ~ holds for all < y < Yq, we have found a Maass cusp form. Typically, 
it is enough to (re-)evaluate for a few different, but well chosen values of y. 

(3) The Laplacian on X is an essentially self-adjoint operator. Hence, all eigenvalues 
X are real. Note, however, that the matrix C(X) is not hermite by construction. As 
such there is no guarantee that the linearisation of C{X + h)a ~ results in real 
solutions. But, if it turns out that a solution becomes independent of y then the 
eigenvalue X = X + h + 0(h 2 ) becomes real, and the Maass form can be made real. 

(4) It may happen that there is no eigenvalue X near a trial value X. In this case, there 
are still 2M(Yq) formal solutions of (C'~ (A)C(A) + hl)a = 0, but \h\ is too large, 
or the formal solutions depend on the height y of the horocycle. 

(5) According to Remark^ a good choice for the height y of the horocycle is cruical. In 
addition, we have to take care that C'(X) is regular and well-conditioned. However, 
it is hard to predict a good value for the height y of the horocycle in advance. 
Moreover, a good choice of y depends on X. 

We fix some reasonable value of y and search with this value of y for Maass 
forms. Sometimes, our choice of y may be good, and sometimes not. Whenever 
the choice of y is not good, we may miss Maass forms. We compensate for this by 
taking another reasonable value of y and run the algorithm again. 

(6) It may happen that eigenvalues degenerate. In principle, this is no problem for 
our algorithm. However, the speed of numerical convergence and the numerically 
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implied errors are strongly affected. Whenever a degenerated eigenvalue occurs, a 
good choice for the height y of the horocycle becomes even more important. 

(7) One may argue that Maass cusp forms, if completely desymmetrised, are conjec- 
tured to be non-degenerated. However, the algorithm does not distinguish whether 
a degeneracy happens between Maass cusp forms, or whether the degeneracy is be- 
tween a Maass cusp form and another formal solution of {C'~ l {\)C{X)-\-hl\a = 0. 
In particular, degeneracies with formal solutions happen frequently. 

(8) In order to speed up the numerics, we do not take the trial values < Ai < A2 < 
... < A max to be sufficiently close together, and hence intentionally risk to oversee 
solutions. We compensate for this by using Weyl's average law and Turing bounds 
after the fact to figure out whether eigenvalues have been overseen. Any missing 
eigenvalues are eventually found with the aid of Algorithm \lfy 

(9) For each trial value X, we search for Maass forms in a neighbourhood. The neigh- 
bourhoods may overlap. Additionally, we may re-run the algorithm. For this reason, 
we typically find each Maass form more than once. Whenever we find a Maass form, 
we consider it only, if it does not match any previously found Maass form. 



6. Turing bound 

From Remark 12 5]) and 12 8]) it is obvious that we may miss some eigenvalues. In the 
current section, we use average Weyl's law and Turing bounds to detect whether eigenvalues 
have been missed, and to figure out small intervals where the missing eigenvalues are to be 
found. 

Let M(t) = #{A I l/4<A<i 2 + 1/4} count the number of eigenvalues which fall into 
the interval [1/4, t 2 + 1/4]. Average Weyl's law reads 

N{t) = M(t) + S(t), 

where M{t) is a smooth approximation to M{t) such that the average value of the Weyl 
remainder S{t) becomes zero in the limit t — > 00, 

(S(t)) :=±jf S(r)dr^?0. 

For Maass forms on certain hyperbolic surfaces, average Weyl's law has been derived 
explicitly [6J, [T7] . An important example is Weyl's average law on the modular surface. 

Theorem 13 (Average Weyl's law on the modular surface [6]). On the modular surface 
X = SL(2, Z)\H we have for t>0: 

. 1 2 2t , t 131 

M(t) = — r — log- 



12" vr"° ev ^72 144 
The proof is given in [6] . Another proof can be found in |17j . 
Remark 14. For X = SL(2,Z)\% there are no small eigenvalues < A < 1/4, |23| [5] . 
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Figure 1. The fluctuations S nnm {t) for SL(2,Z). 



o 6*0 ioo i-:-:-o i moo 

Figure 2. Mean (S num {i)}, with the eigenvalue Ai 359 i 6 = 1/4 + 
1300.0191255 2 removed. 



Let A^ num (i) count the number of numerically found eigenvalues in the interval [1/4, t 2 + 
1/4]. The difference AA num (i) — M(t) is a non-positive integer whose absolute value counts 
the number of solutions which have been overseen numerically. 

As M(t) is unknown, we consider S num (t) = M num (t) - M(t), instead. S num (t) is a 
fluctuating function around its average value, Af nnm (t) — Af(t), see figure [T] From the 
graph of (S num (t)} := j f Q * S num (T)dT, we can read off whether all eigenvalues have been 
found. If all eigenvalues are found, (S nnra (t)) equals (S(t)} and tends to zero in the limit of 
large t. As soon as a solution gets overseen near A = r 2 + l/4 with r > 0, (S num (t)) deviates 
significantly from (S(t)) for t > r, (S num (t)} < (S{t)} + \ £(-l)dr. This is visualised in 
figure[2j where we have intentionally removed the eigenvalue A135916 = 1/4+1300. 0191255 2 . 
From figure [2] we can read off that at least one eigenvalue is missing and that the first 
eigenvalue which is missing is somewhere near A = 1/4 + 1300 2 . 

On the other hand, if we are in doubt whether a given list of eigenvalues is consecutive, 
we can intentionally insert a fake "eigenvalue" near the end of the list. For instance, if 
we intentionally insert a fake "eigenvalue" Af a k e = 1/4 + 1300 2 to the list of consecutive 
eigenvalues, the graph of (S nnia (t\) will clearly show that there is one eigenvalue too much 
near A = 1/4 + 1300 2 , see figure |[ 
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Figure 3. Mean (5 num (t)), with the fake "eigenvalue" A fake = 1/4 + 1300 2 inserted. 

If we know Turing bounds for (S(t)), the comparison of {S nnm (t)) against these bounds 
becomes conclusive. For the modular group, Turing bounds have been derived explicitly 
®- 

Theorem 15 (Turing bounds for the modular group [6]). Consider the modular group 
r = SL(2,Z), and define E(t) = (l + 6.59125/ logt) (vr/(12 logt)) 2 . Then for t > 1, we 
have the lower and upper Turing bounds 

Ei ower {t) := -2E(t) < (S{t)) < E{t) =: E upper {t). 

The proof is given in [6] . 

Theorem 16 (Turing's method |33j). If we add a fake "eigenvalue", \f a ke> near the end 
of a list of eigenvalues, and if (S num (t)) then exceeds the upper Turing bound, the list of 
eigenvalues is consecutive in the interval 1/4 < A < \f a ke- 

The importance of this theorem is highlighted in 21(15]. For the proof, see |33j . 

Theorem 17. Let a list of eigenvalues be given for which we compute {S num (t)). Let lower 
and upper Turing bounds be given such that —1/2 < Ei ower (t) < {S(t)) < E upper (t) < 1/2 
Vt > to, where to > 0. Let 

T := min{r > t | (S num (r)) < E lower (r)} 

and 

t := Tmax (0, (S num (T)) - E upper (T) + l). 

Then, we have: The list of eigenvalues is consecutive for 1/4 < A < t 2 + l/4, but there 
is an eigenvalue missing in t 2 + 1/4 < A<T 2 + l/4. 

Proof. Since there are countably many eigenvalues, the given finite list cannot be complete, 
(S Qum (t)) A/" num (t) - W(t) < -1 < ^iower(t)- T exists and we have (,S num (r)) < 
E\ower(T) < {S(T)), which implies that the list of eigenvalues is not consecutive for 1/4 < 
A<T 2 + l/4. 

If A = r 2 +l/4 is the first missed eigenvalue, we have (5 num (T)) < (S(T)}+± \-l)dr < 
E uppeT (T) + '^. Hence, {S nnm (T)) > E uppcT (T) + t -^ implies t < r and that all eigenvalues 
have been found for 1/4 < A < t 2 + 1/4 if t > 0. □ 

Algorithm 18 (The control algorithm). Let A be a list of numerically found Maass cusp 
forms. Typically, A is empty in the beginning. 
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(1) Let all parameters take reasonable values. 

(2) Determine t and T according to Theorem 17. 

(3) Let t 2 + 1/4 = Ai < A2 < ... < X u = T 2 + 1/4 be trial values, and search for 
eigenvalues in the neighbourhood of each trial value. 

(4) Use Algorithm^ to compute the corresponding Maass cusp forms. 

(5) For each found Maass form, compare whether it is already in the list A. If not, add 
it to the list. 

(6) Let t i d ■= t and T oU := T. 



(7) Determine t and T according to Theorem 17. 

(8) Ift is equal to t id and T is equal to T id, then increase v and change all parameters 
to somewhat different reasonable values. Otherwise, decrease v slightly in order to 
speed up the numerics. 

(9) Continue with^ 



Theoretically, Algorithm 18 runs forever and finds an ever increasing list of consecutive 
Maass cusp forms. 

For groups other than SL(2, Z), Turing bounds have not been derived, yet. It is possible 
to replace Turing bounds by heuristic bounds which follow from intentionally removing 
and inserting eigenvalues near the end of an almost consecutive list of eigenvalues, as was 
demonstrated in figures [2] and |3j However, this requires user interaction and slows down 
the algorithm. Results for certain moonshine groups are published in |17j . 

7. Fluctuations in the Weyl remainder 

Here, we present results for the modular group T = SL(2,Z). It is worth while to 
desymmetrise the Maass cusp forms, as this speeds up their computation. On the modular 
surface X = F\H, Maass cusp forms fall into even, f(—x + iy) = f(x + iy), and odd, 



f{—x + iy) = —f(x + iy), eigenfunctions. Letting Algorithm 18 run for some days, we 
have found 80 443 consecutive even Maass cusp forms with 1/4 < A < 1421. 98 2 + 1/4, but 
missed 1 out of 81 673 with 1/4 < A < 1432. 64 2 + 1/4 by the time when we interrupted the 
algorithm. For the odd symmetry, the algorithm was faster. In the same time, we have 
found 164 232 consecutive odd Maass cusp forms with 1/4 < A < 2000 2 + l/4. The amount 
of data allows us to investigate the asymptotic and statistic properties of the fluctuations 
in the Weyl remainder S(t), numerically. 

From Rudnick |25j we know that the number of eigenvalues obeys a central limit theorem. 
In particular, he considers smoothed windows of length 1/L centered at points t, and 
examins the number Ni,(t) of rj that lie within each such window. In the limit L — > 00 and 
T — > 00 but L = o(logT), the distribution of Ni,(t) as t varies through [T, 2T] approaches 
a Gaussian, 

1 rrr, N L (t)-N£(t) - [* „2 /2 du 

lim - measit <E T, 2T] : — ^ ^ < x} = / e~ u /2 ^=. 

We ask whether the fluctuations in the Weyl remainder (without smoothing) also obey a 
central limit theorem. 
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Figure 4. A graph of ^exp ( - ^ (log log i) 5 / 17 ) \S(t)\ whose magnitude 
of fluctuations slightly increases with t. This eventually lets the limsup of 
the graph diverge in the limit t — > 00. The graph also shows how close the 
lower bound of Li and Sarnak comes to the true asymptotic growth rate. 




Figure 5. A graph of whose magnitude of fluctuations slightly 

decreases with t. This eventually lets the graph vanish in the limit t — > 00. 

In view of figure [IJ we note that the magnitude of the remainder fluctuations grows in t, 
and we may better consider a scaled version of the Weyl remainder. From Li and Sarnak 
[19] we know a lower bound on the asymptotic growth rate of the remainder fluctuations, 

S(t) = fi(^exp(i(loglogt) 5 / 17 )). 

Figure [4] shows how close this lower bound is to the true asymptotic growth rate. Con- 
cerning an upper bound, nothing better than S(t) = 0{t/\ogt) is known analytically. 

Conjecture 19. On the modular surface X = SL(2,Z)\"H, the asymptotic growth rate of 
the remainder fluctuations is bounded from above by 

S{t) = o(Vt). 

Numerical evidence for the conjecture is displayed in figure [5j where the magnitude of 
fluctuations of slightly decreases with t. 

Conjecture 20. Consider the fluctuations in the Weyl remainder S(t) on the modular 
surface X = SL(2, Z)\H. If scaled by log ^ g 1 , the fluctuations obey a central limit theorem. 
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Figure 6. Distribution of l^EMim ^ qM 

curve) as t varies through 
[0, 1421.98], in comparison with the standard Gaussian (dashed curve). 



For any C > 0, we have 

1 log log t ff(t) f x _ u 2 /2 du 

hm — measjt G T, T + L : -= < x\ = / e ' . , 

with a « 0.140. 



A histogram of the scaled distribution is shown in figure 6l The scale factor log ^ g * 

L_J vt 

results from numerical considerations. Namely, if we examine S[t) in the interval [a, b] , the 
mean and the variance of the scaled fluctuations, 

1 f b loglogt 

2 1 y^ loglogt ^ \2 

are numerically independent of i. Quantitative results for /«[t-ioo,t+loo] an d ^ft-ioo.t+ioo] 
in dependence of i are shown in figure [7| We find that the mean vanishes fi ~ and that 
the standard deviation is constant cr ~ 0.140. 

By analogy with the law of the iterated logarithm, one might expect the correct scalings 
for the distribution and extreme values to differ by a factor of y/log log t. However, for the 
S(t) of Riemann zeta, it's thought [9] that they differ by about \/log t, which is quite a 
bit larger. We expect that the difference between the S(t) of Maass form eigenvalues and 
that of Riemann zeros comes from the level statistics. The eigenvalues of the Laplacian 
on arithmetic surfaces are expected to be Poisson distributed [21 13] , whereas the Riemann 
zeros are expected to lie on the critical line and show a distribution in resemblance to the 
eigenvalues of random matrices of the Gaussian unitary ensemble |22l I26j . 

We believe that the iterated logarithm heuristic is closer to the truth. If the Weyl 
remainder would result from a Wiener process, it would strictly follow the law of the 
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Figure 7. Plots of the mean M[t-loo,t+loo] an d the standard deviation 
o"[t-loo,t+ioo] in dependence of t. 



iterated logarithm and we would have 



lim sup ( 



log log t\ 1 / 2 
2t 



\s(t)\ 



a. 



We have checked this numerically and find that the predicted scaling of extreme values 
agrees with our numerical data. The scaling also agrees with the upper bound of Conjecture 
19 However, the numerical values of the extrema are larger than predicted by a factor of 
2.5, i.e. we find 



/loglogfAi/2 
km sup I — - — j \S(t)\ 
t^-oo \ At / 



2.5cj, 



wkich indicates that the fluctuations in the Weyl remainder do not exactly follow the law 
of the iterated lagarithm. 
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